Predation interactions co-occurrence modelling

Packages

Code
# Data manipulation package
library(tidyverse) # Easily Install and Load the 'Tidyverse'
library(psych) # Procedures for Psychological, Psychometric, and Personality Research

# Modelling packages
library(nlme) # Linear and Nonlinear Mixed Effects Models
library(betareg) # Beta Regression
library(glmmTMB) # Generalized Linear Mixed Models using Template Model Builder
library(AICcmodavg) # Model Selection and Multimodel Inference Based on (Q)AIC(c)
library(MuMIn) # Multi-Model Inference

#Modelling utility packages
library(broom) # Convert Statistical Objects into Tidy Tibbles
library(broom.mixed) # Tidying Methods for Mixed Models
library(insight) # Easy Access to Model Information for Various Model Objects
library(modelbased) # Estimation of Model-Based Predictions, Contrasts and Means
library(parameters) # Processing of Model Parameters
library(performance) # Assessment of Regression Models Performance
library(DHARMa) # Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models
library(ggeffects) # Create Tidy Data Frames of Marginal Effects for 'ggplot' from Model Outputs

# Grapich packages
library(GGally) # Extension to 'ggplot2'
library(patchwork) # The Composer of Plots


# Table output packages
library(kableExtra) # Construct Complex Table with 'kable' and Pipe Syntax
library(DT)


library(showtext)
font_add_google("Roboto", "Roboto")
showtext_auto()
showtext_opts(dpi = 300)

Variables description

Table 1: List of variables used to modelling spatial co-occurrence patterns of predator prey interactions
Variable Description Units Source
SIF Spatial species interaction factor. SIF<1: segregation, SIF=1: independence, SIF>1= aggregation index Each study
Ov_coeff Temporal overlap coefficient. Ov_coeff =1: complete activity overlap, Ov_Coeff =0: No activity overlap index Each study
CBMratio Competitors body mass ratio: ln( dominant competitor body mass/ subordinate competitor body mass) index COMBINE data base (Soria et al. 2021)
Dmass_cat Dominant competitor body mass category: <2 kg small, >2 kg and < 5kg Small-Medium, >5 kg and < 15kg Medium-Large, > 15 kg Large category COMBINE data base (Soria et al. 2021)
D_diet Dominant competitor diet category: Carnivore if > 80% of vertebrate diet composition, Herbivore if > 80% of plant diet composition, Insectivore if >80% of invertebrate diet composition, Omnivore rest of the species category PHYLACINE trait data base (Faurby et al. 2020)
Diet_dist Dominant and subordinate competitor category diets: Same if pair of competitors shares de diet category, different if not category PHYLACINE trait data base (Faurby et al. 2020)
Dist_tax Dominant and subordinate competitor taxonomic family distance: Same if the pair of spaces share the taxonomic family, Different if not category COMBINE data base (Soria et al. 2021)
Phy_dist Pairwise mitochondrial DNA phylogenetic distances mitochondrial DNA phylogenetic distances Hassanin et al. (2021)
Abs_lat Absolute latitude of each study when reported coordinates Each study
Continent Continent where the study was carried out category Each study
Label Study labeling category Each study

Data

The data correspond to studies evaluating the spatial and temporal co-occurrence of mammals of the order Carnivora. In this database, mammals of the order Carnivora were considered as predators and mammals of other orders as potential prey.

Spatial co-occurrence data of predators and prey mammals were measured using the species interaction factor (SIF). This is a parameter derived from multi-species occupancy models (Waddle et al. 2010; Richmond, Hines, and Beissinger 2010).

Code
# Load the spatial co-ocurrence data base of predator-prey interactions
spatP_db <- read_delim("Data/Model data/spatP_db.csv", delim = ",")[,-1] %>%
  # Select the variables
  select(SIF, Pred_sp, Pred_family, Pred_mass, Prey_sp, Prey_family, Prey_mass, mass_ratio, P_hunt, Lon, Lat, Lat_abs, Avg_dist, Locality, Label, Pred_evidence) %>% 
  # filter species pair without predator-prey relation
  filter(Pred_evidence == "Yes")
DT::datatable(spatP_db)

The Table 2 contain the general description of the numeric variables of spatial data base. The table was constructed with psych package (Revelle 2022)

Code
Spatnum_summary <- spatP_db %>% 
  select_if(is.numeric) %>% 
  describe(. ,fast = T) 

kbl(Spatnum_summary, caption = "General description of the spatial co-occurrence data base for predator-prey interactions", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
Table 2: General description of the spatial co-occurrence data base for predator-prey interactions.
vars n mean sd min max range se
SIF 1 25 1.081 0.823 0.473 4.656 4.183 0.165
Pred_mass 2 25 37320.851 63890.577 800.000 240500.000 239700.000 12778.115
Prey_mass 3 25 35982.358 46833.523 66.900 135000.000 134933.100 9366.705
mass_ratio 4 25 0.991 2.113 -1.777 5.774 7.551 0.423
Lon 5 23 19.578 91.034 -146.668 113.333 260.002 18.982
Lat 6 23 10.838 36.833 -41.234 63.824 105.057 7.680
Lat_abs 7 23 32.751 18.922 15.402 63.824 48.422 3.946
Avg_dist 8 21 1.214 1.393 0.500 4.000 3.500 0.304

We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient

Code
(SpatP_cor <- select_if(spatP_db, is.numeric) %>% # select numeric variables
  ggpairs(.,
          # Correlation coefficient upper part
        upper = list(continuous= wrap("cor", method= "pearson", 
                                      digits=2, corSize= 80)),
        lower = list( continuous= "smooth")) +
  theme_bw()+
  theme(text=element_text(family = "Roboto")))

Figure 1: Spatial predator-prey co-occurrence corplot

The temporal data corresponds to the co-occurrence of Carnivora mammals with their mammal prey measured by overlap coefficient of kernel activity curves (Ridout and Linkie 2009).

Code
tempP_db <- read_delim("Data/Model data/tempP_db.csv", delim=",")[,-1] %>%
  select(Ov_coeff, Pred_sp, Pred_family, Pred_mass, Prey_sp, Prey_family, Prey_mass, mass_ratio, P_hunt, Lon, Lat, Lat_abs, Samp_dur, Locality, Label,Pred_evidence) %>% 
  # filter species pair without predator-prey relation
  filter(Pred_evidence == "Yes")
DT::datatable(tempP_db )

The Table 3 contain the general description of the numeric variables of spatial data base.

Code
Tempnum_summary <- tempP_db %>% 
  select_if(is.numeric) %>% 
  describe(. ,fast = T) 
  
kbl(Tempnum_summary, caption = "General description of the temporal co-occurrence data base for predator-prey interactions", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
Table 3: General description of the temporal co-occurrence data base for predator-prey interactions.
vars n mean sd min max range se
Ov_coeff 1 693 0.602 0.190 0.040 0.930 0.890 0.007
Pred_mass 2 685 44417.382 41252.700 110.330 240500.000 240389.670 1576.185
Prey_mass 3 693 71389.240 190227.527 36.000 3220000.000 3219964.000 7226.146
mass_ratio 4 685 0.734 1.669 -3.961 5.350 9.311 0.064
Lon 5 693 11.131 79.008 -123.084 149.183 272.267 3.001
Lat 6 693 16.199 21.837 -41.234 56.509 97.743 0.830
Lat_abs 7 693 22.753 14.873 1.686 56.509 54.823 0.565
Samp_dur 8 541 21.182 34.275 1.000 200.000 199.000 1.474

We use the GGAlly package to explore the visual relationship between numerical variables (Schloerke et al. 2021). We used pearson’s correlation coefficient

Code
(TempP_cor <- select_if(tempP_db, is.numeric)  %>% 
  ggpairs(.,
        upper = list(continuous= wrap("cor", method= "pearson", 
                                      digits=2, corSize= 80)),
        lower = list( continuous= "smooth")) +
  theme_bw()+
  theme(text=element_text(family = "Roboto")))

Figure 2: Temporal predator-prey co-occurrence correlogram

Modelling procedure

Standardize numeric covariates

The numerical covariates were standardized to mean 0 and standard deviation 1, which facilitates model fitting and direct comparison of regression coefficients. In addition, pairs of species that do not contain information on any of the variables (cells with NAs) were eliminated.

Because we aimed to evaluate whether the patterns found varied depending on the predator family, we made subsets of the databases. For spatial co-occurrence no subset were made. For temporal co-occurrence we subset the families Felidae and Canidae.

Code
scale_this <- function(data){scale(data) %>% as.vector() }

SP_vars <- c("mass_ratio", "Lat_abs", "Avg_dist")
TP_vars <- c("mass_ratio", "Lat_abs", "Samp_dur")

# All data Spatial
spatP_modsdf <- spatP_db %>%
  drop_na() %>% 
  mutate(across(all_of(SP_vars), scale_this)) 
dim(spatP_modsdf)
[1] 21 16
Code
# All data Temporal
tempP_modsdf <- tempP_db %>%
  drop_na() %>% 
  mutate(across(all_of(TP_vars), scale_this))
dim(tempP_modsdf)
[1] 538  16
Code
# Felidae Temporal
tempP_modsdf_F <- tempP_modsdf %>%
   filter(Pred_family == "Felidae")
dim(tempP_modsdf_F)
[1] 407  16
Code
# Canidae Temporal
tempP_modsdf_C <- tempP_modsdf %>%
   filter(Pred_family == "Canidae")
dim(tempP_modsdf_C)
[1] 85 16

Influential observations

Identify the presence of outliers or observations that may influence the models. For the spatial co-occurrence data we fit a general linear model with all covariates and interactions. We then checked for the presence of extreme data using Laverage distance plots from performance package (Lüdecke et al. 2021).

Code
# Fit model with all data
SPmods_outl <- lm(SIF ~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist, data= spatP_modsdf)

# Check outliers for each model
SP_out <- check_outliers(SPmods_outl)

# Leverage plot for all data model
plot(SP_out)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Figure 3: Spatial co-occurrence leverage distance

Collinearity

We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).

To calculate the VIF we use the Performance package (Lüdecke et al. 2021).

Code
spatP_coll <- lm(SIF ~ mass_ratio+ P_hunt+ Avg_dist, data= spatP_modsdf)

spatP_colltab <- check_collinearity(spatP_coll)


plot(spatP_colltab)+
  labs(subtitle = "All data Spatial model")+
  theme(legend.position = "none")

Figure 4: Spatial co-occurrence collinearity

Model assumptions

We first check the assumptions of the model. To do this we fit the most complex model containing all available variable interactions and random factors. We then visually evaluated the normality of the residuals, linearity and homogeneity of variance using the performance package (Lüdecke et al. 2021).

Code
Spat_diag_m <-lm(SIF ~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist, 
                 random = ~1|Label/Locality,
                 method = "REML",
                 data= spatP_modsdf)

Spat_diagnostic <- check_model(Spat_diag_m, 
                               check = c("homogeneity", "linearity", "qq", "normality", "reqq"))
Spat_diagnostic 

No significant deviations from assumptions were detected.

Random structure

Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).

Code
# Fit model without random effects
Spat_r0 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist,
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)

# Fit model with label as random effect
Spat_r1 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist +(1|Label),
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)

# Fit model with label and Locality as random effects
Spat_r2 <- glmmTMB(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+  Avg_dist +(1|Label/Locality),
                   data=spatP_modsdf, 
                   family=gaussian(link = "identity"), 
                   REML = T)


# Rank the models with the AIC
spat_AICtab <- aictab(cand.set = list(Spat_r0, Spat_r1, Spat_r2),
                       modnames = c("no random",
                                    "Label random", 
                                    "Label/Locality random"),
                                    second.ord = F)

kbl(spat_AICtab, caption = "Random-effects structure selection table for spatial co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
Table 4: Random-effects structure selection table for spatial co-occurrence data
Modnames K AIC Delta_AIC ModelLik AICWt LL Cum.Wt
no random 7 19.478 0 1.000 0.665 -2.739 0.665
Label random 8 21.478 2 0.368 0.245 -2.739 0.910
Label/Locality random 9 23.478 4 0.135 0.090 -2.739 1.000

The Table 4 suggests that models without random effects is best supported. To visualize the variation of the random parameters of each group, we use the estimate_grouplevel function of the modelbased package (Makowski et al. 2020).

Code
estimate_grouplevel(Spat_r1) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

There is no significant variation from each group level in random effects.

Fixed predictor selections

To select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC <2 are considered equally plausible (Burnham and Anderson 2002). Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) (Alain F. Zuur et al. 2009b). because of the amount of data in the spatial co-occurrence model, we only consider in this case additive interactions of the variables. To do this we used the dredge function of the MuMIn package (Barton 2020).

All data

Code
SP_global <- lm(formula = SIF~ mass_ratio+ I(mass_ratio^2)+ P_hunt+ Avg_dist,
                   data=spatP_modsdf, 
                   na.action = "na.fail")


SP_selec <- dredge(SP_global, rank = "AIC", m.lim= c(NA,3))
Code
kbl(SP_selec, 
    caption = "Fixed-effects structure selection table for spatial co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
Fixed-effects structure selection table for spatial co-occurrence data
(Intercept) Avg_dist mass_ratio I(mass_ratio^2) P_hunt df logLik AIC delta weight
15 0.733 NA -0.174 0.082 + 6 8.791 -5.582 0.000 0.27396412
7 0.840 NA -0.098 0.069 NA 4 6.013 -4.026 1.556 0.12582834
1 0.906 NA NA NA NA 2 3.745 -3.490 2.092 0.09626930
11 0.809 NA -0.095 NA + 5 6.614 -3.227 2.355 0.08440608
3 0.906 NA -0.053 NA NA 3 4.466 -2.932 2.650 0.07282897
5 0.878 NA NA 0.029 NA 3 4.073 -2.147 3.435 0.04918139
8 0.839 -0.015 -0.104 0.070 NA 5 6.069 -2.138 3.444 0.04896721
9 0.843 NA NA NA + 4 5.025 -2.049 3.533 0.04683626
12 0.793 0.040 -0.082 NA + 6 6.966 -1.931 3.651 0.04415665
10 0.809 0.063 NA NA + 5 5.853 -1.706 3.876 0.03945156
2 0.906 0.014 NA NA NA 3 3.792 -1.585 3.997 0.03713302
4 0.906 -0.007 -0.056 NA NA 4 4.477 -0.954 4.628 0.02708790
6 0.876 0.018 NA 0.031 NA 4 4.157 -0.315 5.267 0.01967567
13 0.834 NA NA 0.014 + 5 5.106 -0.213 5.369 0.01869998
14 0.803 0.062 NA 0.012 + 6 5.920 0.161 5.743 0.01551355

Confidence intervals explorations

Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).

Code
# Function to get 85%IC
get_ci <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high, Component) %>% 
    filter(Component == "conditional") %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}


get_ci_lm <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high) %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}

# Function to plot

ci_plot <- function(ci_df){
ggplot(ci_df, aes(x=Coefficient, y= Parameter))+
  geom_pointrange(aes(xmin=CI_low, xmax= CI_high,
                      col= Model, linetype= Informative),
                  position = position_dodge2(0.5), linewidth= 1)+
  scale_linetype_manual(values = c("no"="dashed", "yes"= "solid"))+
  geom_vline(xintercept = 0, linetype= "dashed")+
  scale_color_viridis_d()+
  labs(caption = "estimates with 85% ci intervals",
       col= "Model ID")+
  theme_bw()+
    theme(text=element_text(family = "Roboto"))
  
}
Code
# get best models for all data 
SP_best_mods <- get.models(SP_selec, subset = delta <2)

# Apply the function to obtain the table of coefficients of selected the models
SP_best_ci <- map2_df(names(SP_best_mods), SP_best_mods, get_ci_lm) 

ci_plot(SP_best_ci)+ labs(title= "All data")

All data best model 85% CI

Influential observations

Identify the presence of outliers or observations that may influence the models. For the temporal co-ocurrence data we fit a beta family generalized linear model with all covariates and interactions (Cribari-Neto and Zeileis 2010). We then checked for the presence of extreme observations using Cook’s distance plots from performance package (Lüdecke et al. 2021).

Code
# Models for each data
TPmods_outl <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf)

TPmods_outl_F <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf_F)

TPmods_outl_C <-betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, data= tempP_modsdf_C)


TP_out <- check_outliers(TPmods_outl)
TP_out_F <- check_outliers(TPmods_outl_F)
TP_out_C <- check_outliers(TPmods_outl_C)

# Cook's distance plots
plot(TP_out)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

plot(TP_out_F)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

plot(TP_out_C)+
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Outliers of all data temporal model

Outliers of Felidae as predator Temporal model

Outliers of Canidae as predator Temporal model

Temporal co-occurrence cook distance distance

A possible outlier was detected for the Canidae data as predator. It corresponds to observation 27 which exceeded the 0.85 Cook’s distance threshold.

Collinearity

We checked the collinearity of the variables used using the variance inflation factor (VIF). We consider that an VIF value > 10 indicates a high collinearity of the variables (Alain F. Zuur, Ieno, and Elphick 2010).

To calculate the VIF we use the Performance package (Lüdecke et al. 2021).

All species

Code
# Fit models

tempP_coll <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf)
tempP_coll_F <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf_F)
tempP_coll_C <- betareg(Ov_coeff ~mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur, data= tempP_modsdf_C)

# Check collinearity
tempP_colltab <- check_collinearity(tempP_coll)
tempP_colltab_F <- check_collinearity(tempP_coll_F)
tempP_colltab_C <- check_collinearity(tempP_coll_C)

# Vif plots
plot(tempP_colltab)+
  labs(subtitle = "All data Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

plot(tempP_colltab_F)+
  labs(subtitle = "Felidae Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

plot(tempP_colltab_C)+
  labs(subtitle = "Canidae Temporal model")+
  theme_bw()+
  theme(legend.position = "none",
        text=element_text(family = "Roboto"))

(a) VIF of all data model

(b) VIF of Felidae model

(c) VIF of Canidae model

Figure 5: Outliers model detection temporal models

Model assumptions

We verify the model assumptions by visual inspection of the simulated residuals from the DHARMa package (Hartig 2021). To do so we fit the more complex model which includes fixed variables and their interactions, as well as random effects.

All data

Code
Temp_diag <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_diag_res <- simulateResiduals(Temp_diag)
plot(Temp_diag_res)

Since deviations from homogeneity of variance were detected, we evaluated which variables are involved in the model inadequacies.

Code
testCategorical(Temp_diag, tempP_modsdf$P_hunt)
testQuantiles(Temp_diag, tempP_modsdf$mass_ratio)
testQuantiles(Temp_diag, tempP_modsdf$Lat_abs)
testQuantiles(Temp_diag, tempP_modsdf$Samp_dur)
$uniformity
$uniformity$details
catPred: ambush

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.078041, p-value = 0.01688
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: group

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.125, p-value = 0.4413
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: pursuit

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.090531, p-value = 0.398
alternative hypothesis: two-sided


$uniformity$p.value
[1] 0.0168787 0.4413066 0.3979914

$uniformity$p.value.cor
[1] 0.0506361 0.7959827 0.7959827


$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group   2  2.9148 0.05508 .
      535                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.01137
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.01306
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.0075
alternative hypothesis: both

Uniformity of variance test of Hunting strategy

Uniformity of variance test of mass ratio

Uniformity of variance test of absolute latitude

Uniformity of variance test of sampling duration

Code
Temp_diag2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf,
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), 
                   REML = T)

Temp_diag_res2 <- simulateResiduals(Temp_diag2)
plot(Temp_diag_res2)

Felidae

Code
Temp_diag_F <- glmmTMB(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), 
                   REML = T)


Temp_diag_res_F <- simulateResiduals(Temp_diag_F)
plot(Temp_diag_res_F)

Canidae

Code
Temp_diag_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   data=tempP_modsdf_C, 
                   family=beta_family(), 
                   REML = T)

Temp_diag_res_C <- simulateResiduals(Temp_diag_C)
plot(Temp_diag_res_C)

Code
Temp_diag_C2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), 
                   REML = T)

Temp_diag_res_C2 <- simulateResiduals(Temp_diag_C2)
plot(Temp_diag_res_C2)

Code
testCategorical(Temp_diag_C2, tempP_modsdf_C$P_hunt[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$mass_ratio[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$Lat_abs[-27])
testQuantiles(Temp_diag_C2, tempP_modsdf_C$Samp_dur[-27])
$uniformity
$uniformity$details
catPred: group

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.16372, p-value = 0.4185
alternative hypothesis: two-sided

------------------------------------------------------------ 
catPred: pursuit

    One-sample Kolmogorov-Smirnov test

data:  dd[x, ]
D = 0.098545, p-value = 0.6595
alternative hypothesis: two-sided


$uniformity$p.value
[1] 0.4185158 0.6594822

$uniformity$p.value.cor
[1] 0.8370315 0.8370315


$homogeneity
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.5443 0.4628
      82               

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.397
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.9715
alternative hypothesis: both

    Test for location of quantiles via qgam

data:  simulationOutput
p-value = 0.5324
alternative hypothesis: both

Uniformity of variance test of Hunting strategy

Uniformity of variance test of mass ratio

Uniformity of variance test of absolute latitude

Uniformity of variance test of sampling duration

Random structure

Following the Alain F. Zuur et al. (2009a) protocol, we evaluated the relevance of the inclusion of random effects. For this we fit the previously selected model with and without the random effects (Label and Locality of each study). We selected the model with the best structure using the Akaike information criterion (AIC) with the AICcmodavg package (Mazerolle 2023).

Code
# All data
Temp_r0 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2,
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_r1 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label),
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

Temp_r2 <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality),
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   data=tempP_modsdf, 
                   family=beta_family(), REML = T)

temp_AICtab <- aictab(cand.set = list(Temp_r0, Temp_r1, Temp_r2), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random")) 
Code
# Felidae
Temp_r0_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

Temp_r1_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

Temp_r2_F <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

temp_AICtab_F <- aictab(cand.set = list(Temp_r0_F, Temp_r1_F, Temp_r2_F), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random")) 
Code
# Canidae
Temp_r0_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2, 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

Temp_r1_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ Lat_abs + (1|Label), 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

Temp_r2_C <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label/Locality), 
                   data=tempP_modsdf_C[-27,], 
                   family=beta_family(), REML = T)

temp_AICtab_C <- AICcmodavg::aictab(cand.set = list(Temp_r0_C, Temp_r1_C, Temp_r2_C), 
                                  sort = T,
                                  modnames = c("no random",
                                               "Label random",
                                               "Label/Locality random"))
Code
tempP_AICtab_all <- rbind(temp_AICtab, temp_AICtab_F, temp_AICtab_C)

kbl(tempP_AICtab_all, 
    caption = "Random-effects structure selection table for Temporal co-occurrence data", digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F) %>%
  pack_rows("All data",1,3) %>%
  pack_rows("Felidae",4,6) %>% 
  pack_rows("Canidae",7, 9)
Random-effects structure selection table for Temporal co-occurrence data
Modnames K AICc Delta_AICc ModelLik AICcWt LL Cum.Wt
All data
2 Label random 26 -233.714 0.000 1.000 0.752 144.231 0.752
3 Label/Locality random 27 -231.497 2.217 0.330 0.248 144.231 1.000
1 no random 25 -216.630 17.085 0.000 0.000 134.584 1.000
Felidae
21 Label random 18 -252.375 0.000 1.000 0.750 145.069 0.750
31 Label/Locality random 19 -250.174 2.201 0.333 0.250 145.069 1.000
11 no random 17 -232.710 19.665 0.000 0.000 134.141 1.000
Canidae
12 no random 17 19.557 0.000 1.000 0.738 11.858 0.738
22 Label random 18 21.966 2.409 0.300 0.221 12.278 0.959
32 Label/Locality random 19 25.318 5.761 0.056 0.041 12.278 1.000
Code
estimate_grouplevel(Temp_r1) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

estimate_grouplevel(Temp_r1_F) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

estimate_grouplevel(Temp_r1_C) %>% # Get random parameters
  plot()+  #plot
  # aesthetic changes
  theme_bw()+
  theme(text=element_text(family = "Roboto"))

Group level estimates for all temporal data

Group level estimates for Felidae as predator temporal data

Group level estimates for Canidae as predator temporal data

Fixed predictor selections

To select the fixed effects structure we also used the Akaike information criterion (AIC). Models with a delta AIC <2 are considered equally plausible (Burnham and Anderson 2002). Because we are going to choose the fixed effects structure, we refit the previously identified model without using the restricted maximum likelihood (REML) (Alain F. Zuur et al. 2009b). Since co-occurrence patterns may be generated by the interaction of variables, we fit all possible combinations of variables, limiting to a maximum of three parameters per model. To do this we will used the dredge function of the MuMIn package (Barton 2020).

All data

Code
TP_global <- glmmTMB(formula = Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2+ (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), 
                   REML = F,
                    na.action = "na.fail")

TP_selec <- dredge(TP_global, 
                   rank = "AIC", 
                   fixed = c("disp(mass_ratio)","disp(Lat_abs)","disp(Samp_dur)"), 
                   m.lim= c(NA,6))
Code
kbl(TP_selec, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)
Fixed-effects structure selection table for temporal co-occurrence data
cond((Int)) disp((Int)) cond(Lat_abs) cond(mass_ratio) cond(I(mass_ratio^2)) cond(P_hunt) cond(Samp_dur) cond(Lat_abs:mass_ratio) cond(Lat_abs:P_hunt) cond(Lat_abs:Samp_dur) cond(mass_ratio:P_hunt) cond(mass_ratio:Samp_dur) cond(I(mass_ratio^2):Lat_abs) cond(I(mass_ratio^2):mass_ratio) cond(I(mass_ratio^2):P_hunt) cond(I(mass_ratio^2):Samp_dur) cond(P_hunt:Samp_dur) disp(Lat_abs) disp(mass_ratio) disp(Samp_dur) df logLik AIC delta weight
74 0.451 1.911 0.085 NA NA + NA NA + NA NA NA NA NA NA NA NA 0.026 -0.016 -0.071 11 169.361 -316.723 0.000 2.076382e-01
16409 0.461 1.921 NA NA NA + 0.148 NA NA NA NA NA NA NA NA NA + 0.027 -0.002 -0.036 11 169.228 -316.455 0.268 1.816312e-01
26 0.465 1.911 0.083 NA NA + 0.132 NA NA NA NA NA NA NA NA NA NA 0.037 -0.006 -0.047 10 167.737 -315.474 1.248 1.112248e-01
25 0.470 1.910 NA NA NA + 0.137 NA NA NA NA NA NA NA NA NA NA 0.021 0.002 -0.060 9 166.710 -315.421 1.302 1.082803e-01
27 0.471 1.912 NA -0.037 NA + 0.139 NA NA NA NA NA NA NA NA NA NA 0.026 -0.002 -0.059 10 167.205 -314.409 2.314 6.529908e-02
267 0.441 1.921 NA -0.092 NA + NA NA NA NA + NA NA NA NA NA NA -0.007 -0.015 -0.092 11 168.180 -314.359 2.364 6.368705e-02
10 0.442 1.908 0.088 NA NA + NA NA NA NA NA NA NA NA NA NA NA 0.026 -0.006 -0.063 9 165.971 -313.941 2.782 5.167866e-02
9 0.445 1.906 NA NA NA + NA NA NA NA NA NA NA NA NA NA NA 0.006 0.002 -0.079 8 164.873 -313.746 2.977 4.686686e-02
29 0.477 1.911 NA NA -0.009 + 0.137 NA NA NA NA NA NA NA NA NA NA 0.020 0.000 -0.059 10 166.769 -313.539 3.184 4.225497e-02
12 0.442 1.911 0.085 -0.032 NA + NA NA NA NA NA NA NA NA NA NA NA 0.029 -0.009 -0.063 10 166.347 -312.695 4.028 2.770819e-02
11 0.446 1.908 NA -0.035 NA + NA NA NA NA NA NA NA NA NA NA NA 0.010 -0.001 -0.079 9 165.319 -312.637 4.085 2.692625e-02
14 0.450 1.909 0.088 NA -0.009 + NA NA NA NA NA NA NA NA NA NA NA 0.024 -0.008 -0.061 10 166.037 -312.073 4.650 2.030532e-02
13 0.452 1.907 NA NA -0.009 + NA NA NA NA NA NA NA NA NA NA NA 0.004 0.000 -0.078 9 164.928 -311.857 4.866 1.822389e-02
15 0.454 1.909 NA -0.036 -0.010 + NA NA NA NA NA NA NA NA NA NA NA 0.008 -0.005 -0.078 10 165.385 -310.771 5.952 1.058802e-02
4109 0.437 1.910 NA NA 0.004 + NA NA NA NA NA NA NA NA + NA NA -0.002 -0.010 -0.084 11 166.135 -310.270 6.453 8.242532e-03
17 0.384 1.891 NA NA NA NA 0.140 NA NA NA NA NA NA NA NA NA NA 0.030 0.034 -0.065 7 160.425 -306.851 9.872 1.491318e-03
2055 0.370 1.910 NA -0.135 -0.015 NA NA NA NA NA NA NA NA 0.04 NA NA NA 0.036 0.025 -0.091 9 162.046 -306.092 10.631 1.020527e-03
21 0.403 1.894 NA NA -0.020 NA 0.141 NA NA NA NA NA NA NA NA NA NA 0.026 0.027 -0.063 8 160.717 -305.435 11.288 7.347469e-04
18 0.375 1.892 0.043 NA NA NA 0.137 NA NA NA NA NA NA NA NA NA NA 0.037 0.031 -0.061 8 160.710 -305.419 11.304 7.290206e-04
1 0.358 1.888 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0.014 0.034 -0.086 6 158.613 -305.226 11.497 6.618472e-04
19 0.383 1.892 NA -0.017 NA NA 0.141 NA NA NA NA NA NA NA NA NA NA 0.032 0.033 -0.065 8 160.537 -305.074 11.649 6.134844e-04
8213 0.399 1.896 NA NA -0.015 NA 0.167 NA NA NA NA NA NA NA NA -0.03 NA 0.026 0.031 -0.062 9 161.498 -304.997 11.726 5.902105e-04
146 0.367 1.891 0.043 NA NA NA 0.097 NA NA 0.069 NA NA NA NA NA NA NA 0.038 0.032 -0.059 9 161.398 -304.797 11.926 5.340025e-04
22 0.394 1.894 0.045 NA -0.021 NA 0.138 NA NA NA NA NA NA NA NA NA NA 0.033 0.023 -0.059 9 161.023 -304.046 12.676 3.669839e-04
2 0.349 1.890 0.049 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 0.024 0.031 -0.080 7 158.954 -303.907 12.815 3.423555e-04
5 0.377 1.890 NA NA -0.020 NA NA NA NA NA NA NA NA NA NA NA NA 0.011 0.027 -0.083 7 158.897 -303.795 12.928 3.236104e-04
1030 0.377 1.901 0.007 NA -0.032 NA NA NA NA NA NA NA 0.044 NA NA NA NA 0.028 0.022 -0.073 9 160.862 -303.724 12.999 3.123284e-04
23 0.402 1.895 NA -0.018 -0.021 NA 0.142 NA NA NA NA NA NA NA NA NA NA 0.028 0.025 -0.063 9 160.834 -303.667 13.056 3.035919e-04
20 0.374 1.893 0.042 -0.016 NA NA 0.138 NA NA NA NA NA NA NA NA NA NA 0.039 0.031 -0.061 9 160.800 -303.600 13.123 2.935664e-04
3 0.357 1.889 NA -0.016 NA NA NA NA NA NA NA NA NA NA NA NA NA 0.017 0.033 -0.086 7 158.707 -303.413 13.310 2.673875e-04
531 0.382 1.893 NA -0.018 NA NA 0.139 NA NA NA NA 0.009 NA NA NA NA NA 0.032 0.032 -0.065 9 160.568 -303.137 13.586 2.328955e-04
6 0.368 1.893 0.051 NA -0.021 NA NA NA NA NA NA NA NA NA NA NA NA 0.020 0.023 -0.077 8 159.264 -302.527 14.196 1.717009e-04
4 0.348 1.891 0.047 -0.014 NA NA NA NA NA NA NA NA NA NA NA NA NA 0.026 0.030 -0.080 8 159.027 -302.055 14.668 1.355573e-04
7 0.376 1.892 NA -0.017 -0.020 NA NA NA NA NA NA NA NA NA NA NA NA 0.013 0.025 -0.083 8 158.996 -301.992 14.730 1.314138e-04
36 0.342 1.891 0.054 -0.005 NA NA NA -0.047 NA NA NA NA NA NA NA NA NA 0.026 0.023 -0.078 9 159.903 -301.807 14.916 1.197700e-04
8 0.368 1.893 0.049 -0.015 -0.021 NA NA NA NA NA NA NA NA NA NA NA NA 0.022 0.021 -0.077 9 159.342 -300.683 16.039 6.829379e-05

Felidae

Code
TP_Felidae <- glmmTMB(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+ P_hunt+ Lat_abs+ Samp_dur)^2 + (1|Label), 
                   data=tempP_modsdf_F,
                   family=beta_family(),
                   REML = F,
                   na.action = "na.fail")

TP_selec_F <- dredge(TP_Felidae, rank = "AIC", m.lim= c(NA,3))
Code
kbl(TP_selec_F, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data When Felidae is predator", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)
Fixed-effects structure selection table for temporal co-occurrence data When Felidae is predator
cond((Int)) disp((Int)) cond(Lat_abs) cond(mass_ratio) cond(I(mass_ratio^2)) cond(P_hunt) cond(Samp_dur) cond(Lat_abs:mass_ratio) cond(Lat_abs:P_hunt) cond(Lat_abs:Samp_dur) cond(mass_ratio:P_hunt) cond(mass_ratio:Samp_dur) cond(I(mass_ratio^2):Lat_abs) cond(I(mass_ratio^2):mass_ratio) cond(I(mass_ratio^2):P_hunt) cond(I(mass_ratio^2):Samp_dur) cond(P_hunt:Samp_dur) df logLik AIC delta weight
18 0.464 + 0.132 NA NA NA 0.105 NA NA NA NA NA NA NA NA NA NA 5 155.852 -301.703 0.000 0.072760138
26 0.478 + 0.122 NA NA + 0.098 NA NA NA NA NA NA NA NA NA NA 6 156.839 -301.677 0.026 0.071818780
10 0.466 + 0.127 NA NA + NA NA NA NA NA NA NA NA NA NA NA 5 155.749 -301.499 0.204 0.065698485
2 0.451 + 0.139 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 4 154.711 -301.422 0.281 0.063221822
20 0.464 + 0.123 -0.053 NA NA 0.110 NA NA NA NA NA NA NA NA NA NA 6 156.662 -301.323 0.380 0.060168047
12 0.466 + 0.119 -0.050 NA + NA NA NA NA NA NA NA NA NA NA NA 6 156.495 -300.991 0.713 0.050951490
4 0.450 + 0.131 -0.050 NA NA NA NA NA NA NA NA NA NA NA NA NA 5 155.434 -300.869 0.835 0.047935630
267 0.461 + NA -0.076 NA + NA NA NA NA + NA NA NA NA NA NA 6 156.253 -300.506 1.197 0.039984181
8213 0.448 + NA NA 0.011 NA 0.150 NA NA NA NA NA NA NA NA -0.055 NA 6 156.203 -300.406 1.297 0.038040265
36 0.443 + 0.142 -0.060 NA NA NA -0.058 NA NA NA NA NA NA NA NA NA 6 156.172 -300.344 1.359 0.036875703
27 0.474 + NA -0.062 NA + 0.108 NA NA NA NA NA NA NA NA NA NA 6 155.997 -299.993 1.710 0.030943502
22 0.477 + 0.132 NA -0.014 NA 0.104 NA NA NA NA NA NA NA NA NA NA 6 155.941 -299.881 1.822 0.029257739
74 0.468 + 0.123 NA NA + NA NA + NA NA NA NA NA NA NA NA 6 155.874 -299.747 1.956 0.027361915
146 0.466 + 0.132 NA NA NA 0.110 NA NA -0.012 NA NA NA NA NA NA NA 6 155.872 -299.744 1.959 0.027321172
25 0.473 + NA NA NA + 0.102 NA NA NA NA NA NA NA NA NA NA 5 154.844 -299.688 2.016 0.026559184
19 0.460 + NA -0.062 NA NA 0.118 NA NA NA NA NA NA NA NA NA NA 5 154.843 -299.685 2.018 0.026527940
14 0.478 + 0.127 NA -0.013 + NA NA NA NA NA NA NA NA NA NA NA 6 155.830 -299.660 2.043 0.026196281
6 0.465 + 0.139 NA -0.016 NA NA NA NA NA NA NA NA NA NA NA NA 5 154.829 -299.657 2.046 0.026161193
17 0.460 + NA NA NA NA 0.112 NA NA NA NA NA NA NA NA NA NA 4 153.731 -299.461 2.242 0.023719992
11 0.461 + NA -0.059 NA + NA NA NA NA NA NA NA NA NA NA NA 5 154.701 -299.403 2.300 0.023034874
9 0.461 + NA NA NA + NA NA NA NA NA NA NA NA NA NA NA 4 153.671 -299.341 2.362 0.022338752
8 0.463 + 0.132 -0.049 -0.014 NA NA NA NA NA NA NA NA NA NA NA NA 6 155.527 -299.055 2.648 0.019354848
3 0.445 + NA -0.059 NA NA NA NA NA NA NA NA NA NA NA NA NA 4 153.449 -298.897 2.806 0.017890622
1 0.446 + NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 3 152.445 -298.891 2.812 0.017831724
1030 0.457 + 0.108 NA -0.007 NA NA NA NA NA NA NA 0.034 NA NA NA NA 6 155.232 -298.464 3.239 0.014402868
16409 0.476 + NA NA NA + 0.101 NA NA NA NA NA NA NA NA NA + 6 155.102 -298.204 3.499 0.012649437
29 0.485 + NA NA -0.013 + 0.102 NA NA NA NA NA NA NA NA NA NA 6 154.916 -297.833 3.871 0.010505669
23 0.471 + NA -0.060 -0.013 NA 0.117 NA NA NA NA NA NA NA NA NA NA 6 154.914 -297.829 3.874 0.010485735
531 0.459 + NA -0.063 NA NA 0.114 NA NA NA NA 0.014 NA NA NA NA NA 6 154.910 -297.821 3.883 0.010442640
21 0.473 + NA NA -0.015 NA 0.112 NA NA NA NA NA NA NA NA NA NA 5 153.837 -297.674 4.030 0.009702304
15 0.472 + NA -0.058 -0.012 + NA NA NA NA NA NA NA NA NA NA NA 6 154.762 -297.523 4.180 0.008999102
13 0.474 + NA NA -0.014 + NA NA NA NA NA NA NA NA NA NA NA 5 153.761 -297.522 4.181 0.008992944
5 0.461 + NA NA -0.017 NA NA NA NA NA NA NA NA NA NA NA NA 4 152.578 -297.157 4.546 0.007492844
7 0.459 + NA -0.058 -0.015 NA NA NA NA NA NA NA NA NA NA NA NA 5 153.547 -297.094 4.609 0.007263183
4109 0.477 + NA NA -0.019 + NA NA NA NA NA NA NA NA + NA NA 6 154.047 -296.094 5.609 0.004404904
2055 0.458 + NA -0.048 -0.014 NA NA NA NA NA NA NA NA -0.004 NA NA NA 6 153.559 -295.118 6.585 0.002704093

Canidae

Code
tempP_modsdf_C2 <- slice(tempP_modsdf_C, -27)

TP_Canidae <- betareg(Ov_coeff ~ (mass_ratio+ I(mass_ratio^2)+P_hunt + Lat_abs+Samp_dur)^2,
                         data=tempP_modsdf_C2, 
                         na.action = "na.fail")

TP_selec_C <- dredge(TP_Canidae, rank = "AIC", m.lim= c(NA,3))
Code
kbl(TP_selec_C, 
    caption = "Fixed-effects structure selection table for temporal co-occurrence data When Canidae is predator", 
    digits = 3) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), 
                full_width = F)
Fixed-effects structure selection table for temporal co-occurrence data When Canidae is predator
(Intercept) Lat_abs mass_ratio I(mass_ratio^2) P_hunt Samp_dur Lat_abs:mass_ratio Lat_abs:P_hunt Lat_abs:Samp_dur mass_ratio:P_hunt mass_ratio:Samp_dur I(mass_ratio^2):Lat_abs I(mass_ratio^2):mass_ratio I(mass_ratio^2):P_hunt I(mass_ratio^2):Samp_dur P_hunt:Samp_dur df logLik AIC delta weight
74 -0.336 0.539 NA NA + NA NA + NA NA NA NA NA NA NA NA 5 20.551 -31.102 0.000 0.391806844
36 0.191 0.100 0.221 NA NA NA -0.231 NA NA NA NA NA NA NA NA NA 5 20.088 -30.176 0.926 0.246623687
1030 0.323 -0.071 NA -0.142 NA NA NA NA NA NA NA 0.12 NA NA NA NA 5 19.023 -28.047 3.055 0.085067398
26 -0.101 0.180 NA NA + 0.226 NA NA NA NA NA NA NA NA NA NA 5 17.832 -25.663 5.438 0.025833351
18 0.098 0.176 NA NA NA 0.203 NA NA NA NA NA NA NA NA NA NA 4 16.510 -25.019 6.082 0.018721043
22 0.178 0.187 NA -0.053 NA 0.209 NA NA NA NA NA NA NA NA NA NA 5 17.334 -24.667 6.435 0.015697255
16409 0.016 NA NA NA + 0.483 NA NA NA NA NA NA NA NA NA + 5 17.224 -24.447 6.654 0.014064257
25 0.016 NA NA NA + 0.184 NA NA NA NA NA NA NA NA NA NA 4 16.206 -24.413 6.689 0.013822389
2055 0.263 NA -0.232 -0.045 NA NA NA NA NA NA NA NA 0.061 NA NA NA 5 17.138 -24.276 6.826 0.012909008
1 0.186 NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA 2 14.044 -24.088 7.013 0.011753387
17 0.207 NA NA NA NA 0.157 NA NA NA NA NA NA NA NA NA NA 3 14.997 -23.994 7.107 0.011213252
20 0.123 0.163 0.072 NA NA 0.227 NA NA NA NA NA NA NA NA NA NA 5 16.970 -23.940 7.162 0.010911143
9 0.015 NA NA NA + NA NA NA NA NA NA NA NA NA NA NA 3 14.960 -23.921 7.181 0.010808779
2 0.100 0.133 NA NA NA NA NA NA NA NA NA NA NA NA NA NA 3 14.960 -23.920 7.182 0.010803512
10 -0.071 0.135 NA NA + NA NA NA NA NA NA NA NA NA NA NA 4 15.905 -23.810 7.291 0.010228099
146 0.128 0.162 NA NA NA 0.200 NA NA 0.075 NA NA NA NA NA NA NA 5 16.841 -23.682 7.420 0.009592431
19 0.227 NA 0.087 NA NA 0.193 NA NA NA NA NA NA NA NA NA NA 4 15.669 -23.339 7.763 0.008079372
6 0.175 0.143 NA -0.050 NA NA NA NA NA NA NA NA NA NA NA NA 4 15.668 -23.336 7.766 0.008068637
5 0.258 NA NA -0.045 NA NA NA NA NA NA NA NA NA NA NA NA 3 14.600 -23.199 7.902 0.007534955
21 0.282 NA NA -0.046 NA 0.160 NA NA NA NA NA NA NA NA NA NA 4 15.594 -23.188 7.914 0.007492804
29 0.092 NA NA -0.036 + 0.184 NA NA NA NA NA NA NA NA NA NA 5 16.569 -23.138 7.963 0.007309855
14 0.013 0.142 NA -0.040 + NA NA NA NA NA NA NA NA NA NA NA 5 16.374 -22.748 8.353 0.006013944
3 0.196 NA 0.056 NA NA NA NA NA NA NA NA NA NA NA NA NA 3 14.331 -22.662 8.439 0.005761653
13 0.091 NA NA -0.036 + NA NA NA NA NA NA NA NA NA NA NA 4 15.318 -22.637 8.465 0.005688479
27 0.053 NA 0.033 NA + 0.193 NA NA NA NA NA NA NA NA NA NA 5 16.276 -22.552 8.550 0.005451713
4109 0.227 NA NA -0.099 + NA NA NA NA NA NA NA NA + NA NA 5 16.264 -22.529 8.573 0.005388523
23 0.289 NA 0.077 -0.039 NA 0.192 NA NA NA NA NA NA NA NA NA NA 5 16.110 -22.220 8.881 0.004619106
4 0.113 0.124 0.040 NA NA NA NA NA NA NA NA NA NA NA NA NA 4 15.103 -22.206 8.895 0.004586076
11 0.016 NA 0.001 NA + NA NA NA NA NA NA NA NA NA NA NA 4 14.961 -21.921 9.181 0.003976818
12 -0.105 0.141 -0.026 NA + NA NA NA NA NA NA NA NA NA NA NA 5 15.948 -21.896 9.205 0.003927985
7 0.260 NA 0.045 -0.040 NA NA NA NA NA NA NA NA NA NA NA NA 4 14.789 -21.578 9.523 0.003350397
8 0.180 0.137 0.025 -0.047 NA NA NA NA NA NA NA NA NA NA NA NA 5 15.724 -21.448 9.654 0.003138977
531 0.228 NA 0.088 NA NA 0.194 NA NA NA NA 0.002 NA NA NA NA NA 5 15.669 -21.339 9.763 0.002972685
8213 0.292 NA NA -0.052 NA 0.196 NA NA NA NA NA NA NA NA -0.023 NA 5 15.631 -21.263 9.839 0.002861382
15 0.087 NA -0.004 -0.036 + NA NA NA NA NA NA NA NA NA NA NA 5 15.319 -20.639 10.463 0.002094630
267 0.125 NA 0.097 NA + NA NA NA NA + NA NA NA NA NA NA 5 15.182 -20.364 10.737 0.001826176

Confidence intervals explorations

Code
# Function to get 85%IC
get_cibeta <- function(mod_name, mods){
  ci_df <- parameters::model_parameters(mods, ci = 0.85) %>% 
    select(Parameter, Coefficient, CI_low, CI_high) %>% 
    filter(Parameter != "(Intercept)") %>%
    mutate(Model = mod_name,
           Informative = case_when(
    CI_low < 0 & CI_high < 0 ~ "yes",
    CI_low > 0 & CI_high > 0 ~ "yes",
    TRUE ~ "no")) 
  return(ci_df)   
}

Consistent with the model selection strategy (AIC) and to avoid uninformative variables, we explored the 85% confidence intervals of the regression coefficients of the best models (Sutherland et al. 2023). Variables whose 85% CI overlap 0 are considered as uninformative and therefore not included in the inference (Arnold 2010).

Code
# get best models for all data 
TP_best_mods <- get.models(TP_selec, subset = delta <2, REML = T)
# get best models for felidae subset data
TP_best_mods_F <- get.models(TP_selec_F, subset = delta <2, REML = T)
# get best models for Canidae subset data
TP_best_mods_C <- get.models(TP_selec_C, subset = delta <2)


# Apply the function to obtain the table of coefficients of selected the models
TP_best_ci <- map2_df(names(TP_best_mods), TP_best_mods, get_ci) 
TP_best_ci_F <- map2_df(names(TP_best_mods_F), TP_best_mods_F, get_ci)
TP_best_ci_C <- map2_df(names(TP_best_mods_C), TP_best_mods_C, get_cibeta)


ci_plot(TP_best_ci)+ labs(title= "All data")
ci_plot(TP_best_ci_F)+ labs(title= "Felidae as predator")
ci_plot(TP_best_ci_C)+ labs(title= "Canidae as predator")

All data best model 85% CI

Felidae data best model 85% CI

Canidae data best model 85% CI

Summary selected models

Note that some variables have coefficients whose 85% confidence intervals overlap 0. These variables are retained when they have interactions with other variables and these interactions do not overlap with 0.

Spatial co-occurrence models

Code
SP_mod1 <- lm(formula = SIF~ mass_ratio +I(mass_ratio^2)+  P_hunt,
                   data=spatP_modsdf)

SP_mod2 <- lm(formula = SIF~ mass_ratio +I(mass_ratio^2),
                   data=spatP_modsdf)


SP_mods <- list(SP_mod1, SP_mod2) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("SIF~ mass_ratio +I(mass_ratio)^2+P_hunt", 5),
                  rep("SIF~ mass_ratio +I(mass_ratio)^2", 3)),
         Component= "conditional",
         Effects= "fixed") %>% 
  select(-SE,-t,-df_error)

kbl(SP_mods, 
    caption = "Spatial co-occurrence selected models", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F)
Spatial co-occurrence selected models
Parameter Coefficient CI CI_low CI_high p Model Component Effects
(Intercept) 0.733 0.85 0.624 0.843 0.000 SIF~ mass_ratio +I(mass_ratio)^2+P_hunt conditional fixed
mass_ratio -0.174 0.85 -0.276 -0.073 0.020 SIF~ mass_ratio +I(mass_ratio)^2+P_hunt conditional fixed
I(mass_ratio^2) 0.082 0.85 0.018 0.147 0.073 SIF~ mass_ratio +I(mass_ratio)^2+P_hunt conditional fixed
P_huntgroup -0.014 0.85 -0.267 0.239 0.935 SIF~ mass_ratio +I(mass_ratio)^2+P_hunt conditional fixed
P_huntpursuit 0.221 0.85 0.069 0.374 0.043 SIF~ mass_ratio +I(mass_ratio)^2+P_hunt conditional fixed
(Intercept) 0.840 0.85 0.753 0.927 0.000 SIF~ mass_ratio +I(mass_ratio)^2 conditional fixed
mass_ratio -0.098 0.85 -0.174 -0.021 0.072 SIF~ mass_ratio +I(mass_ratio)^2 conditional fixed
I(mass_ratio^2) 0.069 0.85 0.008 0.130 0.108 SIF~ mass_ratio +I(mass_ratio)^2 conditional fixed

Goftest of spatial selected models

Code
SP_best_goft_list <- list(SP_mod1, SP_mod2) 
SP_best_goft <- lapply(SP_best_goft_list, check_model, check = c("homogeneity", "linearity", "qq", "normality", "reqq"))

Temporal co-occurrence models

Code
TP_mod1 <-  glmmTMB(formula = Ov_coeff ~  Lat_abs *P_hunt + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mod2 <-  glmmTMB(formula = Ov_coeff ~ P_hunt* Samp_dur + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mod3 <-  glmmTMB(formula = Ov_coeff ~ P_hunt+ Samp_dur + (1|Label), 
                   data=tempP_modsdf, 
                   dispformula = ~ mass_ratio+ Samp_dur+Lat_abs,
                   family=beta_family(), REML = T)

TP_mods <- list(TP_mod1, TP_mod2, TP_mod3) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~  Lat_abs *P_hunt + (1|Label)", 11),
                  rep("Ov_coeff ~ P_hunt* Samp_dur + (1|Label)", 11),
                  rep("Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)", 9))) %>% 
  select(-SE, -z, -df_error, -Group)
Code
TP_mod1_F <-  glmmTMB(formula = Ov_coeff ~  Lat_abs + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mod2_F <-  glmmTMB(formula = Ov_coeff ~  Lat_abs + Samp_dur+ (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)  

TP_mod3_F <-  glmmTMB(formula = Ov_coeff ~ mass_ratio *P_hunt + (1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T) 

TP_mod4_F <-  glmmTMB(formula = Ov_coeff ~  I(mass_ratio^2)* Samp_dur+(1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mod5_F <-  glmmTMB(formula = Ov_coeff ~  mass_ratio+ Samp_dur +(1|Label), 
                   data=tempP_modsdf_F, 
                   family=beta_family(), REML = T)

TP_mods_F <- list(TP_mod1_F, TP_mod2_F, TP_mod3_F, TP_mod4_F, TP_mod5_F) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~  Lat_abs + (1|Label)", 4),
                  rep("Ov_coeff ~  Lat_abs + Samp_dur+ (1|Label)", 5),
                  rep("Ov_coeff ~ mass_ratio * P_hunt + (1|Label)", 6),
                  rep("Ov_coeff ~  I(mass_ratio^2) * Samp_dur+(1|Label)", 6),
                  rep("Ov_coeff ~  mass_ratio+ Samp_dur +(1|Label)", 5))) %>% 
  select(-SE, -z, -df_error,-Group)  
Code
TP_mod1_C  <-betareg(Ov_coeff ~ P_hunt * Lat_abs,
                         data=tempP_modsdf_C2)

TP_mod2_C  <-betareg(Ov_coeff ~ mass_ratio * Lat_abs,
                         data=tempP_modsdf_C2)

TP_mods_C <- list(TP_mod1_C, TP_mod2_C) %>% 
  map(model_parameters, digits= 2, ci= 0.85) %>% 
  reduce(rbind) %>%
  mutate(Model= c(rep("Ov_coeff ~ P_hunt * Lat_abs", 4),
                  rep("Ov_coeff ~ mass_ratio * Lat_abs", 4)),
         Component= "conditional",
         Effects= "fixed") %>% 
  select(-SE, -z, -df_error)
Code
TP_mods_tab <- rbind(TP_mods, TP_mods_F, TP_mods_C)

kbl(TP_mods_tab, 
    caption = "Temporal co-occurrence selected models", digits = 3) %>%
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = F) %>% 
    pack_rows("Temporal co-occurrence models -all data",1, nrow(TP_mods)) %>%
    pack_rows("Temporal co-occurrence models -Felidae as predator",1+ nrow(TP_mods),nrow(TP_mods)+ nrow(TP_mods_F)) %>%
  pack_rows("Temporal co-occurrence models -Canidae as predator",1+nrow(TP_mods)+ nrow(TP_mods_F), nrow(TP_mods)+ nrow(TP_mods_F)+ nrow(TP_mods_C))
Temporal co-occurrence selected models
Parameter Coefficient CI CI_low CI_high p Effects Component Model
Temporal co-occurrence models -all data
(Intercept) 0.444 0.85 0.354 0.535 0.000 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs 0.083 0.85 -0.012 0.178 0.209 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
P_huntgroup -0.564 0.85 -0.769 -0.360 0.000 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
P_huntpursuit -0.148 0.85 -0.318 0.022 0.209 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs:P_huntgroup 0.249 0.85 0.034 0.463 0.095 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs:P_huntpursuit -0.146 0.85 -0.293 0.001 0.153 fixed conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
(Intercept) 1.911 0.85 1.822 2.000 0.000 fixed dispersion Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
mass_ratio -0.016 0.85 -0.108 0.075 0.798 fixed dispersion Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Samp_dur -0.072 0.85 -0.166 0.021 0.264 fixed dispersion Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
Lat_abs 0.028 0.85 -0.075 0.130 0.695 fixed dispersion Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
SD (Intercept) 0.313 0.85 0.238 0.411 NA random conditional Ov_coeff ~ Lat_abs *P_hunt + (1|Label)
(Intercept) 0.453 0.85 0.358 0.548 0.000 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntgroup -0.438 0.85 -0.642 -0.235 0.002 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntpursuit -0.154 0.85 -0.309 0.001 0.153 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Samp_dur 0.147 0.85 0.037 0.258 0.054 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntgroup:Samp_dur 0.202 0.85 -0.101 0.505 0.338 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
P_huntpursuit:Samp_dur -0.162 0.85 -0.294 -0.030 0.077 fixed conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
(Intercept) 1.919 0.85 1.830 2.009 0.000 fixed dispersion Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
mass_ratio -0.002 0.85 -0.092 0.089 0.975 fixed dispersion Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Samp_dur -0.045 0.85 -0.134 0.045 0.473 fixed dispersion Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
Lat_abs 0.024 0.85 -0.077 0.125 0.731 fixed dispersion Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
SD (Intercept) 0.333 0.85 0.257 0.432 NA random conditional Ov_coeff ~ P_hunt* Samp_dur + (1|Label)
(Intercept) 0.463 0.85 0.368 0.557 0.000 fixed conditional Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
P_huntgroup -0.468 0.85 -0.667 -0.269 0.001 fixed conditional Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
P_huntpursuit -0.169 0.85 -0.324 -0.014 0.116 fixed conditional Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Samp_dur 0.135 0.85 0.030 0.240 0.064 fixed conditional Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
(Intercept) 1.910 0.85 1.822 1.999 0.000 fixed dispersion Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
mass_ratio 0.001 0.85 -0.089 0.091 0.988 fixed dispersion Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Samp_dur -0.062 0.85 -0.152 0.027 0.316 fixed dispersion Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Lat_abs 0.022 0.85 -0.077 0.122 0.745 fixed dispersion Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
SD (Intercept) 0.329 0.85 0.256 0.423 NA random conditional Ov_coeff ~ P_hunt+ Samp_dur + (1|Label)
Temporal co-occurrence models -Felidae as predator
(Intercept) 0.444 0.85 0.354 0.535 0.000 fixed conditional Ov_coeff ~ Lat_abs + (1|Label)
Lat_abs 0.137 0.85 0.042 0.231 0.038 fixed conditional Ov_coeff ~ Lat_abs + (1|Label)
(Intercept) 7.977 0.85 7.205 8.831 NA fixed dispersion Ov_coeff ~ Lat_abs + (1|Label)
SD (Intercept) 0.320 0.85 0.245 0.417 NA random conditional Ov_coeff ~ Lat_abs + (1|Label)
(Intercept) 0.458 0.85 0.368 0.548 0.000 fixed conditional Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
Lat_abs 0.129 0.85 0.036 0.223 0.046 fixed conditional Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
Samp_dur 0.101 0.85 -0.001 0.204 0.154 fixed conditional Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
(Intercept) 7.973 0.85 7.202 8.826 NA fixed dispersion Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
SD (Intercept) 0.309 0.85 0.235 0.408 NA random conditional Ov_coeff ~ Lat_abs + Samp_dur+ (1|Label)
(Intercept) 0.454 0.85 0.365 0.544 0.000 fixed conditional Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
mass_ratio -0.075 0.85 -0.136 -0.015 0.073 fixed conditional Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
P_huntgroup -0.244 0.85 -0.591 0.102 0.310 fixed conditional Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
mass_ratio:P_huntgroup 0.316 0.85 0.056 0.576 0.080 fixed conditional Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
(Intercept) 7.966 0.85 7.188 8.827 NA fixed dispersion Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
SD (Intercept) 0.307 0.85 0.230 0.410 NA random conditional Ov_coeff ~ mass_ratio * P_hunt + (1|Label)
(Intercept) 0.443 0.85 0.342 0.544 0.000 fixed conditional Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
I(mass_ratio^2) 0.010 0.85 -0.040 0.060 0.777 fixed conditional Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
Samp_dur 0.147 0.85 0.042 0.253 0.045 fixed conditional Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
I(mass_ratio^2):Samp_dur -0.053 0.85 -0.089 -0.018 0.032 fixed conditional Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
(Intercept) 7.997 0.85 7.220 8.858 NA fixed dispersion Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
SD (Intercept) 0.316 0.85 0.240 0.416 NA random conditional Ov_coeff ~ I(mass_ratio^2) * Samp_dur+(1|Label)
(Intercept) 0.454 0.85 0.361 0.546 0.000 fixed conditional Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
mass_ratio -0.061 0.85 -0.120 -0.002 0.139 fixed conditional Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
Samp_dur 0.115 0.85 0.012 0.219 0.108 fixed conditional Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
(Intercept) 7.990 0.85 7.215 8.847 NA fixed dispersion Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
SD (Intercept) 0.326 0.85 0.250 0.425 NA random conditional Ov_coeff ~ mass_ratio+ Samp_dur +(1|Label)
Temporal co-occurrence models -Canidae as predator
(Intercept) -0.336 0.85 -0.596 -0.076 0.063 fixed conditional Ov_coeff ~ P_hunt * Lat_abs
P_huntpursuit 0.675 0.85 0.351 0.999 0.003 fixed conditional Ov_coeff ~ P_hunt * Lat_abs
Lat_abs 0.539 0.85 0.305 0.774 0.001 fixed conditional Ov_coeff ~ P_hunt * Lat_abs
P_huntpursuit:Lat_abs -0.620 0.85 -0.909 -0.331 0.002 fixed conditional Ov_coeff ~ P_hunt * Lat_abs
(Intercept) 0.191 0.85 0.029 0.353 0.089 fixed conditional Ov_coeff ~ mass_ratio * Lat_abs
mass_ratio 0.221 0.85 0.088 0.354 0.017 fixed conditional Ov_coeff ~ mass_ratio * Lat_abs
Lat_abs 0.100 0.85 -0.040 0.240 0.305 fixed conditional Ov_coeff ~ mass_ratio * Lat_abs
mass_ratio:Lat_abs -0.231 0.85 -0.337 -0.124 0.002 fixed conditional Ov_coeff ~ mass_ratio * Lat_abs

Goftest of temporal selected models

Code
TP_best_goft_list <- list(TP_mod1, TP_mod2, TP_mod3) 
TP_best_goft <- lapply(TP_best_goft_list, simulateResiduals, plot=T)

Ov_coeff ~ Lat_abs * P_hunt + (1 | Label) goftest -all data

Ov_coeff ~ P_hunt * Samp_dur + (1 | Label) goftest - -all data

Ov_coeff ~ P_hunt + Samp_dur + (1 | Label) goftest- -all data

Code
TP_best_goft_list_F <- list(TP_mod1_F, TP_mod2_F, TP_mod3_F, TP_mod4_F, TP_mod5_F) 
TP_best_goft_F <- lapply(TP_best_goft_list_F, simulateResiduals, plot=T)

Ov_coeff ~ Lat_abs + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ Lat_abs + Samp_dur + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ mass_ratio * P_hunt + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ I(mass_ratio^2) * Samp_dur + (1 | Label) goftest -Felidae as predator

Ov_coeff ~ mass_ratio + Samp_dur + (1 | Label) goftest -Felidae as predator

Code
TP_best_goft_list_C <- list(TP_mod1_C, TP_mod2_C) 
TP_best_goft_C <- lapply(TP_best_goft_list_C, plot)

Predictions

Code
# Spatial data

SP_pred_mass <- ggemmeans(SP_mod1, terms = c("mass_ratio[all]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(spatP_db$mass_ratio), "scaled:scale"))+ attr(scale(spatP_db$mass_ratio), "scaled:center"))

SP_pred_hunt <- ggemmeans(SP_mod1, terms = c("P_hunt"),
                          ci_level=.85,
                          ci.lvl = .85)

# Temporal data

TP_pred_lat1 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.10:2.35 by=.2]",
                                    "P_hunt[pursuit]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat2 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.3:2.36 by=.2]",
                                    "P_hunt[ambush]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat3 <- ggemmeans(TP_mod1, 
                          terms = c("Lat_abs[-1.15:2.37 by=.2]",
                                    "P_hunt[group]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lat <- rbind(TP_pred_lat1, TP_pred_lat2, TP_pred_lat3)

# PhuntxSampling
TP_pred_sxph1 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.58:2.57 by=.2]",
                                    "P_hunt[pursuit]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph2 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.58:5.20 by=.2]",
                                    "P_hunt[ambush]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph3 <- ggemmeans(TP_mod2, 
                          terms = c("Samp_dur[-0.53:2.14 by=.2]",
                                    "P_hunt[group]"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

TP_pred_sxph <- rbind(TP_pred_sxph1, TP_pred_sxph2, TP_pred_sxph3)

# Felidae temporal data
#laitud
TP_pred_lat_F <- ggemmeans(TP_mod1_F, terms = c("Lat_abs"),
                          ci_level=.85,
                          ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

# masa * P_hunt

TP_pred_mxhunt_F1 <- ggemmeans(TP_mod3_F, 
                               terms = c("mass_ratio[-2.41:2.56 by=.2]", 
                                         "P_hunt[ambush]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"))

TP_pred_mxhunt_F2 <- ggemmeans(TP_mod3_F, 
                               terms = c("mass_ratio[-1.48:1.57 by=.2]", 
                                         "P_hunt[group]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"))


TP_pred_mxhunt_F <- rbind(TP_pred_mxhunt_F1, TP_pred_mxhunt_F2)

# Sampdur

TP_pred_samp_F <- ggemmeans(TP_mod5_F, 
                               terms = c("Samp_dur[-0.58:5.20 by=0.2]"),
                               ci_level=.85,
                               ci.lvl = .85) %>% 
  mutate(samp_real= (x* attr(scale(tempP_db$Samp_dur), "scaled:scale"))+ attr(scale(tempP_db$Samp_dur), "scaled:center"))

# Canidae temporal data

TP_pred_lxh_C1 <-  ggemmeans(TP_mod1_C, 
                            terms = c("Lat_abs[-1.1:1.7 by=.2]", 
                                      "P_hunt[pursuit]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lxh_C2 <-  ggemmeans(TP_mod1_C, 
                            terms = c("Lat_abs[-1.1:2.3 by=.2]", 
                                      "P_hunt[group]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(lat_real= (x* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center"))

TP_pred_lxh_C <- rbind(TP_pred_lxh_C1, TP_pred_lxh_C2)


TPlat_vec_C <- c(max(tempP_modsdf_C2$Lat_abs), mean(tempP_modsdf_C2$Lat_abs),
               min(tempP_modsdf_C2$Lat_abs))

TP_pred_mxl_C <-  ggemmeans(TP_mod2_C, terms = 
                              c("mass_ratio", 
                                "Lat_abs[TPlat_vec_C]"),
                            ci_level=.85,
                            ci.lvl = .85) %>% 
  mutate(mass_real= (x* attr(scale(tempP_db$mass_ratio), "scaled:scale"))+ attr(scale(tempP_db$mass_ratio), "scaled:center"),
         lat_real= (as.numeric(as.character(group))* attr(scale(tempP_db$Lat_abs), "scaled:scale"))+ attr(scale(tempP_db$Lat_abs), "scaled:center")) %>% 
  mutate(across(lat_real, round, 2)) 

Plots

Code
(SP_mass_plot <- ggplot(SP_pred_mass)+
    geom_hline(yintercept = 1, linetype= "dashed", size= 1)+
    geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , alpha=0.7, fill= "#fde725")+
  geom_line(aes(x= mass_real, y=predicted),
            linewidth= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Spatial Overlap",
       tag= "A")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code
(SP_phunt_plot <- ggplot(SP_pred_hunt)+
   geom_hline(yintercept = 1, linetype= "dashed", size= 1)+
   geom_errorbar(aes(x= x, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , size= 1, width= 0.1)+
  geom_point(aes(x= x, y=predicted),
            size= 4)+
  labs(x= "Hunt strategy",
        y= "Spatial Overlap",
       tag= "B")+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"))
  )

Code
# Temporal models

(TP_latxp_plot <- ggplot(TP_pred_lat)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group,
                  group= group)
              , alpha= 0.7, linewidth= 1)+
    scale_fill_viridis_d()+
  geom_line(aes(x= lat_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
       group= "Hunt strategy",
       linetype= "Hunt strategy",
       fill= "Hunt strategy",
       tag= "C")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="vertical",
         legend.position = c(0.7, 0.20),
         legend.background = element_blank()))

Code
(TP_sampxp_plot <- ggplot(TP_pred_sxph)+
   geom_ribbon(aes(x= samp_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group,
                  group= group)
              , alpha= 0.7, linewidth= 1)+
    scale_fill_viridis_d()+
  geom_line(aes(x= samp_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Survey duration (months)",
        y= "Temporal Overlap",
       group= "Hunt strategy",
       linetype= "Hunt strategy",
       fill= "Hunt strategy",
       tag= "H")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="vertical",
         legend.position = c(0.8, 0.20),
         legend.background = element_blank()))

Code
# Felidae latitude
(TP_lat_plot_F <- ggplot(TP_pred_lat_F)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high),
               alpha=0.7, fill= "#fde725")+
   geom_line(aes(x= lat_real, y=predicted),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
        tag= "E")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code
(TP_mxh_plot_F <- ggplot(TP_pred_mxhunt_F)+
   geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group, 
                  group= group,
                  linetype= group)
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= mass_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Temporal Overlap",
       group= "Hunt \nstrategy",
       fill= "Hunt \nstrategy",
       linetype= "Hunt \nstrategy",
       tag= "D")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code
(TP_samp_plot_F <- ggplot(TP_pred_samp_F)+
    geom_ribbon(aes(x= samp_real, y= predicted,
                  ymin= conf.low, ymax= conf.high)
              , alpha=0.7, fill= "#fde725")+
  geom_line(aes(x= samp_real, y=predicted),
            size= 0.8 )+
  labs(x= "Survey duration (months)",
        y= "Temporal Overlap",
       tag= "I")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto")))

Code
(TP_lxh_plot_C <- ggplot(TP_pred_lxh_C)+
   geom_ribbon(aes(x= lat_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= group, 
                  group= group,
                  linetype= group)
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= lat_real, y=predicted, group= group,
                linetype= group),
            size= 0.8 )+
  labs(x= "Absolute latitude",
        y= "Temporal Overlap",
       group= "Hunt \nstrategy",
       fill= "Hunt \nstrategy",
       linetype= "Hunt \nstrategy",
       tag= "F")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction="horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code
(TP_mxl_plot_C <- ggplot(TP_pred_mxl_C)+
   geom_ribbon(aes(x= mass_real, y= predicted,
                  ymin= conf.low, ymax= conf.high,
                  fill= as.factor(lat_real), 
                  group= as.factor(lat_real),
                  linetype= as.factor(lat_real))
              , alpha=0.7)+
    scale_fill_viridis_d()+
  geom_line(aes(x= mass_real, y=predicted, group= as.factor(lat_real),
                linetype= as.factor(lat_real)),
            size= 0.8 )+
  labs(x= "ln(Mass ratio)",
        y= "Temporal Overlap",
       group= "Absolute \nlatitude",
       fill= "Absolute \nlatitude",
       linetype= "Absolute \nlatitude",
       tag= "G")+
  scale_x_continuous(expand = c(0,0))+
  theme_bw()+
   theme(text = element_text(size=13, family = "Roboto"),
         legend.direction= "horizontal",
         legend.position = c(0.55, 0.10),
         legend.background = element_blank()))

Code
(PPrediction_plot <- (SP_mass_plot+ SP_phunt_plot+ TP_latxp_plot+ 
                     TP_mxh_plot_F+ TP_lat_plot_F+ TP_lxh_plot_C+ 
                       TP_mxl_plot_C+ TP_sampxp_plot+ TP_samp_plot_F)+
  plot_layout(ncol = 3))

Code
ggsave(filename = "Figs/Ppreds_plot.png", plot = PPrediction_plot,
        width = 12, height = 10)

References

Arnold, Todd W. 2010. “Uninformative Parameters and Model Selection Using Akaike’s Information Criterion.” The Journal of Wildlife Management 74 (6): 1175–78. https://doi.org/10.2193/2009-367.
Barton, Kamil. 2020. “MuMIn: Multi-Model Inference.” https://CRAN.R-project.org/package=MuMIn.
Burnham, Kenneth P., and David R. Anderson. 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd ed. New York: Springer-Verlag. //www.springer.com/us/book/9780387953649.
Cribari-Neto, Francisco, and Achim Zeileis. 2010. “Beta Regression in r 34. https://doi.org/10.18637/jss.v034.i02.
Faurby, Søren, Rasmus Ø Pedersen, Matt Davis, Simon D. Schowanek, Scott Jarvie, Alexandre Antonelli, and Jens-Christian Svenning. 2020. MegaPast2Future/PHYLACINE_1.2: PHYLACINE Version 1.2.1. Zenodo. https://doi.org/10.5281/zenodo.3690867.
Hartig, Florian. 2021. “DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models.” https://CRAN.R-project.org/package=DHARMa.
Hassanin, Alexandre, Géraldine Veron, Anne Ropiquet, Bettine Jansen van Vuuren, Alexis Lécu, Steven M. Goodman, Jibran Haider, and Trung Thanh Nguyen. 2021. “Evolutionary History of Carnivora (Mammalia, Laurasiatheria) Inferred from Mitochondrial Genomes.” PLOS ONE 16 (2): e0240770. https://doi.org/10.1371/journal.pone.0240770.
Lüdecke, Daniel, Mattan S. Ben-Shachar, Indrajeet Patil, Philip Waggoner, and Dominique Makowski. 2021. Performance: An r Package for Assessment, Comparison and Testing of Statistical Models” 6: 3139. https://doi.org/10.21105/joss.03139.
Makowski, Dominique, Mattan S. Ben-Shachar, Indrajeet Patil, and Daniel Lüdecke. 2020. “Estimation of Model-Based Predictions, Contrasts and Means.” https://github.com/easystats/modelbased.
Mazerolle, Marc J. 2023. “AICcmodavg: Model Selection and Multimodel Inference Based on (q)AIC(c).” https://cran.r-project.org/package=AICcmodavg.
Revelle, William. 2022. “Psych: Procedures for Psychological, Psychometric, and Personality Research.” https://CRAN.R-project.org/package=psych.
Richmond, Orien M. W., James E. Hines, and Steven R. Beissinger. 2010. “Two-Species Occupancy Models: A New Parameterization Applied to Co-Occurrence of Secretive Rails.” Ecological Applications 20 (7): 2036–46. https://doi.org/10.1890/09-0470.1.
Ridout, M. S., and M. Linkie. 2009. “Estimating Overlap of Daily Activity Patterns from Camera Trap Data.” Journal of Agricultural, Biological, and Environmental Statistics 14 (3): 322–37. https://doi.org/10.1198/jabes.2009.08038.
Schloerke, Barret, Di Cook, Joseph Larmarange, Francois Briatte, Moritz Marbach, Edwin Thoen, Amos Elberg, and Jason Crowley. 2021. “GGally: Extension to ’Ggplot2’.” https://CRAN.R-project.org/package=GGally.
Soria, Carmen D., Michela Pacifici, Moreno Di Marco, Sarah M. Stephen, and Carlo Rondinini. 2021. “COMBINE: A Coalesced Mammal Database of Intrinsic and Extrinsic Traits.” Ecology 102 (6): e03344. https://doi.org/10.1002/ecy.3344.
Sutherland, C., D. Hare, P. J. Johnson, D. W. Linden, R. A. Montgomery, and E. Droge. 2023. “Practical Advice on Variable Selection and Reporting Using AIC.” Proceedings of the Royal Society B: Biological Sciences. https://ora.ox.ac.uk/objects/uuid:38650aad-20a6-4980-af4f-7ded0968e908?fbclid=IwAR29QmTXxKFknl9BI_SKw-6I5gL-9NmdaytggzFomgMOPMY9V50q9ohJQ6k.
Waddle, J.Hardin, Robert M. Dorazio, Susan C. Walls, Kenneth G. Rice, Jeff Beauchamp, Melinda J. Shuman, and Frank Mazzotti. 2010. “A New Parameterization for Estimating Co-Occurrence of Interacting Species.” Ecological Applications 20 (5): 302315. https://doi.org/https://doi.org/10.1890/09-0850.1.
Zuur, Alain F., Elena N. Ieno, and Chris S. Elphick. 2010. “A Protocol for Data Exploration to Avoid Common Statistical Problems.” Methods in Ecology and Evolution 1 (1): 3–14. https://doi.org/https://doi.org/10.1111/j.2041-210X.2009.00001.x.
Zuur, Alain F, Elena N Ieno, Neil J Walker, Anatoly A Saveliev, and Graham M Smith. 2009a. Mixed Effects Models and Extensions in Ecology with r. Vol. 574. Springer.
———. 2009b. Mixed Effects Models and Extensions in Ecology with r. 1st ed. Statistics for Biology and Health. Springer.